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Abstract 

The confinement of critical fluctuations in soft media induces critical Casimir forces acting on the 
confining surfaces. The temperature and geometry dependences of such forces are characterized 
by universal scaling functions. A novel approach is presented to determine them for films via 
Monte Carlo simulations of lattice models. The method is based on an integration scheme of free 
energy differences. Our results for the Ising and the XY universality class compare favourably 
with corresponding experimental results for wetting layers of classical binary liquid mixtures and 
of ^He, respectively. 
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I. INTRODUCTION 



Recent progress in understanding the features of effective forces induced by confined 
fluctuations, botli quantum and thermal, reveals the potential relevance of these so-called 
Casimir forces [l| for numerous applications, ranging from microelectromechanical systems 
(MEMS) to the physics of colloids j^, ^, 4]. Thermal fluctuation-induced Casimir forces fc 
acting on the confining surfaces of fluids near critical points jsl are of particular interest 
because they become largely independent of the microscopic details of the system, acquiring 
a universal character [sl,!^,!?], and they can be switched on and off upon varying, e.g., the 
temperature. Moreover, by changing the surface chemistry they can be relatively easily 
turned from attractive to repulsive ^, in contrast to the Casimir force stemming from 
electromagnetic fluctuations, for which such a possibility is currently debated as being very 
desirable to avoid stiction in MEMS, but difficult to achieve. Finite-size scaling theory 
(see, e.g., ref. jsl) predicts that the temperature dependence of the critical Casimir force 
fc is described by universal scaling functions which depend on the bulk universality class 
(UC) of the confined medium and on the surface UCs of the confining surfaces j^. The 
latter are related to the boundary conditions (BC) ^ imposed by the surfaces on the 
relevant fluctuating field, i.e., the order parameter (OP) of the underlying second-order phase 
transition. In spite of intensive theoretical and experimental efforts, the current knowledge 
of these scaling functions is still rather limited even for relevant UCs such as the Ising one, 
which characterizes the critical behaviour of simple fluids and binary liquid mixtures. In 
three spatial dimensions (3D) the only available results refer, theoretically, to films with 
periodic BC (PBC), investigated via Monte Carlo (MC) simulations Q], or field-theoretical 
methods (Dirichlet, Neumann BC, PBC) [3], and, experimentally, to complete wetting films 
of binary liquid mixture belonging to the surface UC characterized by symmetry-breaking 
surface fields jo]. The corresponding BC (H — ) of opposing surface fields reflect the fact that 
the two confining surfaces exhibit opposite adsorption preferences for the two species of the 



mixture. At the bulk critical point the dependence of fc on the thickness L of the 



ilm 



turned out to be in good agreement with the corresponding theoretical predictions [10|, 
However, the determination of the full temperature dependence of fc from these very difficult 
experiments suffers from significant statistical and systematic uncertainties, enhancing the 
need of theoretical insight. Indeed, several features of the associated scaling function, such 
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as its global shape and its dependence on the spatial dimensionality d and BC still await 
theoretical investigations. Exact results are available in d = 2 [l^ and d > 4 [l^ (mean- 
field theory) both for (H — ) and (++) BC, the latter corresponding to the case in which 
both confining surfaces exhibit preference for the same species of the mixture. Proposals 
to measure the temperature dependence of the Casimir force between a colloid and a fiat 
surface or between colloidal particles dissolved in a near-critical binary liquid mixture ^ 
call for a detailed theoretical analysis of the associated scaling behaviour. The relevant 
missing pieces, mentioned above, in the theoretical analysis of the scaling behaviour of fc 
require to account for the fiuctuations, including the dimensional crossover occurring in a 
film. This is a rather challenging task, especially if the OP profile is inhomogeneous across 
the film, as in the cases we are interested in. With these elements out of reach of current 
analytical techniques, MC simulations provide a useful alternative approach. The available 
MC results for the d = 3 Ising model are restricted to the case of PBC 8|, in which the 
scaling function of fc can be determined — up to a normalization factor — by numerically 
measuring the expectation value of a suitable lattice stress tensor. The purpose of the 
present contribution is to present a novel approach for the MC simulation of the Casimir 
force and to provide data for the scaling behaviour of fc- We focus on the Ising UC with 
the experimentally relevant BC (++) and (H — ). We also compare our results with those in 
ref. js], providing an independent test of the method proposed therein. Our method is based 
on an integration scheme of free energy differences and it has the advantage, compared to 
the latter, of providing the absolute value for fc and of being applicable for arbitrary BC. 
The comparison with the experimental data in ref. reveals good agreement. 

Measurements [l^ of the equilibrium thickness of ^He wetting films near the superfiuid 
temperature Tx provide an experimental determination of the scaling function of fc for the 
XY UC with Dirichlet BC on both surfaces, corresponding to the so-called ordinary surface 
UC gI. These BC are due to the fact that the superfiuid OP vanishes at the surfaces. In 
this case more analytical and numerical results are available. For temperatures T >Tx field- 
theoretical calculations [3] of the scaling function are in agreement with the experimental 
data jl3], whereas its behaviour for T <^ T\ has been determined by accounting for He- 
specific features related to capillary- wavelike surface fiuctuations [l^. In addition, valuable 
information on the shape of the scaling function in the critical region has been obtained on 
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the basis of Landau- Ginzburg theory [15[. Recent MC simulations [16| have nicely confirmed 
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and extended the available analytic and experimental results. Since our approach differs from 



16| we also present our results for the scaling function in this case, providing 



the one in ref. 
a valuable test. 

In a film geometry with thickness L and large transverse area A, the Casimir force /c 
per unit area A and in units of k^T = (3~^ is defined as fc{f3,L) = —df^^/dL, where 
= (3L[f — f^^^^iP)] is the excess free energy (which depends on the BC), / is the 
free energy of the film per unit volume V = LA and f^^^^ is the bulk free energy density. 



According to finite-size scaling [3[ the Casimir force takes the universal scaling form 

/c(/?,L)=L-'^^(r(L/e^-)V^) (1) 

where the scaling function -d^x) depends on the BC, t = {Pc — P) / P is the reduced temper- 
ature and ^ = ^^\t\~'^ is the bulk correlation length which controls the spatial exponential 
decay of correlations. The critical exponent u equals 0.6301(4) and 0.662(7) for the 3D Ising 
and XY UCs, respectively are nonuniversal amplitudes above (+) and below (— ) 

T 

c- 

II. COMPUTATION OF THE SCALING FUNCTIONS 

We compute the Casimir force for the Ising and XY models defined on a 3D simple 
cubic lattice in slab geometry (L^ x Ly x with = Ly ^ = L and A = x Ly) 
via the Hamiltonian H = —JYl{ij)^i ' ^i' where the sum is taken over all nearest 

neighbour pairs of sites i and j on the lattice. In the Ising model, Sj has only one component 
Si G {+1,-1}, whereas in the XY model Sj is a two-component vector with modulus |sj| = 1. 
With the Hamiltonian H one finds pc = 0.2216544(3) and = 0.501(2) M] for the Ising 
model, whereas pc = 0.45420(2) and = 0.498(2) [IJ], for the XY model. ^ Temperatures 
and energies are measured in units of J and C,q in units of the lattice spacing. In the x and 
y directions we assume PBC whereas in the z direction we consider periodic, free, and fixed 
BC (i.e., for the Ising model, Sj = -|-1 (-I-) or Si = —1 (— ) at the boundaries). For large 
A, the total free energy F{P, L, A) of such systems decomposes as F{P, L, A) = ALf = 
A[Lf^''^^{p) + p-^f''''{p,L)], where t^'iP.L = oo) is the contribution to F due to the two 
isolated surfaces, macroscopically far apart from each other, whereas f^^iP, L) — f^^{P, oo) 
is the L-dependent finite-size contribution we are interested in. On the lattice ("), fc{P,L) 
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is given by 

where AF(/3, L, A) = F{(3, L, A) - F{/3, L-1,A). 

In general MC methods do note lend themselves to the efficient computation of quantities 
such as F, which cannot be expressed as suitable ensemble averages. However, the free 
energy difference AF{j3, L, A) we are interested in can be cast in such a form via the so- 
called "coupling parameter approach" (see, e.g., ref. 20|). This is a viable alternative to 
the method used in ref. 8| in which AF has been expressed as the ensemble average of a 
lattice stress tensor, which so far is only applicable for PBC. We consider two lattice models 
with the same configuration space C but different Hamiltonians Hi and H2, so that their 
free energies are given by Fi = — ^ In exp(— where indicates the sum over all 
possible configurations belonging to C. F2 — F1 can be conveniently computed by introducing 
the crossover Hamiltonian 

H,,{\) = {l-X)Hi + XH2 (3) 

which depends on the coupling parameter A G [0, 1] and interpolates between Hi and H2 as 
A increases from to 1. Accordingly, the free energy -Fcr(A) = — ln^^exp{—j3Hcr{X)) of a 
system with Hamiltonian iJcr(A) and configuration space C interpolates between Fi and F2. 
The derivative of -Fcr('^) with respect to the coupling parameter, 

dFc.(A) _ Ec(-g2-i/i)e-^^"(^) _ 

takes the form of the canonical ensemble average (. . ■)hci-{x) (with Hamiltonian Hcr{X)) of 
AH = H2 — Hi and therefore it can be efficiently computed via MC simulations. A straight- 
forward integration over A yields the expression for the free energy difference 

F2-Fi= / dA (Aif)^^^(,) = / (5) 
Jo 

in terms of an ensemble average (see, e.g., ref. [2^). 

The Casimir force is related to the difference AF{(3,L,A) (see eq. ([2])) between the free 
energies F{i3,L,A) and F{i3,L — 1,A) of the same model on two lattices with different 
numbers of sites and therefore different configuration spaces. In order to apply the method 
described above for the computation of AF(/3, L, A) one identifies the initial Hamiltonian 
Hi and the configuration space C with the corresponding ones of the model on the lattice 
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(a) Hi (b) H2 



(c) Hcr{X) 



FIG. 1: Bond arrangement for the computation of the free energy difference in eq. ([5|) (see main 
text). 

A X L, as depicted in fig. [T]^a). Accordingly, Fi{P,L,A) = F{P,L,A). The configuration 
space of the final system can be arranged to be equal to C by adding to the model on the 
lattice A X {L — 1) we are actually interested in a two-dimensional lattice of size A with 
suitable degrees of freedom and lateral PBC (see fig. [11(b)). The final Hamiltonian H2 is 
then constructed such that the added layer does not interact with the remaining part of the 
system and therefore F2(/3, L, A) = F{(3, L-1,A) + F2d(/5, A), where F2Z)(/3, A) is the free 
energy of the isolated two-dimensional layer. Although the argument is quite general, we 
focus now on the case of interest in which the lattice degrees of freedom (e.g., spins of the 
Ising model) interact only with their nearest neighbours on the same lattice, with a coupling 
strength J = 1 (indicated by solid bonds in figs. [1] (a) and (b)). The crossover Hamiltonian 
ifcr(A) (see eq. ([3])) additionally depends on the position ko G {1,2, (along the z- 

direction) of the two-dimensional layer which decouples from the rest of the system upon 
passing from A = to A = 1, i.e., from fig. [T](a) to (b). The resulting Hcr{X) is characterized 
by the coupling constants depicted in fig. [Tt^c) whereas AH (see eq. ([5])) can be determined 
as AH = Hcr{\ = 1) — Hcr{\ = 0). AF (see eqs. ([2]) and ([5])) can be finally expressed as 
AF(/3, L, A) = -I{/3, L, A)+F2Di(3, A) from which one has still to subtract f^''^^i/3) in order 
to determine the Casimir force in a slab of thickness L — 1/2 (see eq. ([2])). However, it is 
numerically more convenient to avoid the computation of f^^^^{(3) by considering, instead, 
the difference between the Casimir forces in slabs of thicknesses Li and L2 > Li\ 

A/c(/5, Li, L2, A) ^ fciP, L, -\.A)- /c(/5, - ^, A) 

= /?A-i[J(/5,Li,A)-J(/5,L2,A)], 

in which the contributions of both /*^""^(/5) and F2£)(/3, A) actually cancel. Below we de- 
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scribe the method used to determine in eq. ([T]) on the basis of the numerical data for 
A/c(/3, Li, L2, A). In passing, we note that ahhough iJcr(A) (see fig. [^c)), AH, and there- 
fore {AH)h^^(^x) depend on the choice of ko, /q^ dA {AH)h^^(^x) is actually independent of it, 
as long as the boundary conditions are not affected by the extraction of the ko-th layer. In 
particular, imposing the BC at the boundary layers in the 2-direction, this requires kg ^ 1,L 
for fixed and open BC, whereas for PBC there is no restriction and indeed translational in- 
variance implies that (AH) h^,{\) is actually independent of /cq- In our simulations we have 
taken k^ = L/2. 

Within the MC simulations we compute the ensemble averages (AH) Hc^{\) for different 
values of /?, lattice sizes, and A. Then, via numerical integration (Simpsons method with 20 
points) we calculate the integral J(/5, L, A) in eq. and thus Afc{P, Li,L2, A) (see eq. (jSD) 
with (Li,L2) = (L, 2L) and L = 13,16,20 for the Ising model, whereas L = 10,15,20 for 
the XY model. For a given pair of thicknesses (Li, L2), fixed A and BC, the scaling function 
of the Casimir force can be extracted from the temperature dependence of A/c by using 
the fact that, for large Li^2 and A, fc in eq. scales according to eq. 1^. In particular, it 
is useful to focus on the quantity 

g{y; L,, L2, A) = {L, - 1/2)" A/c(/5(y; L,), L,, L,, A) , (7) 

as a function of y, where P{y; Li) = Pc/[^ + y(-^i ^ 1/2)"-'^/'^] and d = 3; g is expected to 
scale as (see eq. ([T])) 

g{y; Li, L2, A) = e{y) - a-'e{a^'^y) , (8) 

where a = {L2 — l/2)/(Li — 1/2) is the width ratio and d{y) is the lattice estimate of 
9{y) = "^iy / {^qY^^)- For given g eq. ([HD can be solved for 9{y) via an iterative method. 
In the first step one takes O^iy) = g{y; Li, L2, A) as a first approximation of the actual 
6. In turn, this approximant can be improved by taking into account that eq. ([8]) yields 
6{y) = Oo{y) + a~'^9{Q^/^y) ~ O^iy) + a~''-6o{a^^'^y), so that a better approximant Oi{y) 
is provided by Oi{y) = 9o{y) + a~^9o{a^^'^y) where the value of 'do at the point a^/'^y is 
obtained by cubic spline interpolation of the available data. In the expression for 9i one 
can replace 6q by using eq. ([H]), yielding Oiijj) = 6{y) — a~'^'^9{a'^/'^y), which, in turn, can 
be solved as already done before for eq. (IHl) by introducing 6*2 (y) = ^i(y) + a~'^^9i{a'^^'^y) = 
d{y)- a~^'^e{a^/''y), and so on. This iterative procedure yields a sequence of approximants 
^k>i{y) = ^k-i{y) +0;^^'' ^'^9k-i{ci^'' ^^"y)-, which converges very rapidly because the correc- 
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tion to the k-th. approximant is of the order of a ^ i.e., exponentially small in 2''. Already 
for k = 5 one has a~^'^ ^'^ ~ 3.5 x 10^^''' in 3D with a ~ 2( 24]). Accordingly, we approximate 



9(y) = 9k~^Qoiy) by 95{y). The result for the universal scaling function '&{x) = (^{x{^qY^'^) of 
/c, as obtained from a specific pair of lattices with thicknesses (L, 2L), should be indepen- 
dent of the actual value of L, at least for large L. However, for the thicknesses we used in 
our MC simulations, corrections to the leading scaling behaviour are actually relevant 
and affect both the scaling variable 
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x^T{L/^+Y/%l + g^L~-) (9) 

and the scaling function i}{x) which additionally depends on L~'^' : fc{f3,L,A) = 
L"'^d{x,L~'^') ~ L~'^{}{x)[l + L~^'(f){x) + . . .] for large L. In eq. ([9]), u is the leading 
correction-to-scaling exponent u ^ 0.84(4) and 0.79 171] for the Ising and XY UC, respec- 
tively, in 3D. uj' controls the leading corrections to scaling of the estimator fc and its value, 
typically equal to u, can be increased by using suitably improved Hamiltonians and ob- 



servables 17] so that these corrections are reduced. However, we point out that in the 

Jl 

present case surface operators J6| might even yield u' < u. For small lattice sizes, next-to- 
leading corrections to scaling (e.g., ~ L~'^) might also be of relevance, resulting in effective 
L-dependent exponents. A detailed analysis of all these corrections and the determination of 
4>{x) and uj' is beyond the scope of the present paper and requires the study of much larger 
lattices. As a phenomenological ansatz for the effective corrections we take tu' ~ u; ~ 1 
and Q 

/c(/3, L, A) = L-\l + g,L-')-'^{x) , (10) 

which, for suitable choices of gi and g^j, yields a good data collapse of the curves corre- 
sponding to different sizes. We point out that an equally satisfactory data collapse can 
be obtained — within the range of the scaling variable x explored here — by assuming 
a different functional form for the corrections. The resulting estimate of ■^{x) is slightly 
affected by this choice and only larger scale simulations can provide an estimate of -^Ix) 
which is unbiased in this respect. In addition to these corrections to scaling, the simulation 
data depend on the aspect ratio o = L/^/A. Whereas in the case of the XY model this 



dependence is quite pronounced 16|], for the Ising model with (++) and (H — ) BC and for 
X > —6, Afc{P, L, 2L, A) exhibits only a very weak dependence on p already for p < 1/6, as 
we have tested by considering lattices with 1/p = 6, 10, 14. The results we present here for 

8 



^or the XY model we have accounted 



16| by considering multiphcative cor- 



the Ising model refer to lattices with fixed p = 1/6. 
for corrections due to p 7^ in accordance with ref. 
rections 1 + rifP' and (1 + r2p^)~^ to x and /c, respectively, which allow the extrapolation 
of the data on lattices with 1/p = 4, 6, 8 to p ^ 0. 

The computation of the canonical average (AH) HcrW been carried out via an hybrid 
MC method in which a suitable mixture of Wolff and Metropolis algorithms is used 2l|. In 
particular, for the Ising model each hybrid MC step consists of one flip of a Wolff cluster 
according to the Wolff algorithm, typically followed by 9A attempts to flip a spin Sx,y,2 
with 2; G {ko — 1, ko, ko + 1}, which are accepted according to the Metropolis rate 21|. An 
analogous method, with a suitable implementation of Metropolis and Wolff algorithms has 
been used fo. .Ke XY ™odel Q. We tested our MC program successfully by comparing 
g{y, 5, 10, 9) with corresponding transfer-matrix data. 



III. 3D XY MODEL 



We have determined the scaling function of the Casimir force in the 3D XY model with 
free BC, which is of relevance for the ^He experiment mentioned in the introduction. The 
resulting scaling function is plotted in flg. [2l In order to achieve scaling we have accounted 
for corrections to scaling according to eq. ( |T0|) with Qi = 6.4(2), g^^ = 2.1(2) and for the 
corrections due to p 7^ in accordance with ref. 16| with ri = 2.3(2) and r2 = 1.1(1) 



fsee ref. 



16| for details). The scaling function in flg. [2] is compatible, within the errorbars. 



with the one determined in ref. 



161], providing an independent test both of the results 



16|, fc is 



presented there and of our method to compute it. In the approach used in ref. 
computed via the internal energy density u and an integration over the temperature, whereas 
here this is carried out via the free energy density / and an integration over the coupling 
A. In ref. [l^ one takes advantage of the possibly available numerical knowledge of the 
bulk energy density u^^^^ of the model of interest whereas here the analogous information 
on jg rgqui];g(i for the determination of making our approach applicable also 



to cases in which there is no detailed knowledge of u^'^^^ and f^^^^. T 
■t}{x) in flg. [2] compare also very well with the experimental data in ref. 



le MC results for 



13| (we have used 



the experimental value = 1.432 A for the normalization of x). In particular with 

Xmin = —5.29(7) and ^rnin = ''^(a^min) = —1-41(2), it captures properly the corresponding 
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FIG. 2: Scaling function of the Casimir force for the 3D XY bulk UC and so-called ordinary 
surface UC corresponding to free boundary conditions. Our MC data compare very well with the 
corresponding experimental data from ref. [2] (solid line). 

experimental values x^^^^ = —5.7(5) and 1!^^^^^ = —1.30(3) for the pronounced minimum. 
In comparing ■t^min with ^^^^ one has to take into account that, as pointed out above, 
the numerical determination of is actually influenced by the choice of the ansatz for 
the corrections to scaling. Indeed, replacing the multiplicative correction (1 + giL~^)~^ in 
eq. ( fTOl) with one of the form [1 — giL"^) (which is equivalent to the previous one for large 
L), the resulting scaling function displays a good data collapse for ^1 = 3.1(1). It has the 
same shape as the one in fig. [2] but its overall amplitude is reduced by a factor R ~ 0.89. 
With this caveat, our estimates for Xmin and i^min are compatible also with those of ref. 13] 
(—5.3(1) and —1.35(3), respectively). 



IV. 3D ISING MODEL 



For this bulk UC we have determined the scaling functions for three different BC: (H — ), 
(++) (pairs of lattices with sizes (L, 2L) and L = 13, 16, 20) and PBC {L = 10, 16, 20). The 
results are reported in figs. [31 [Hand 0, respectively. Each data point has been averaged over 
at least 10^ hybrid MC steps. The scaling function in fig. [3]has been obtained accounting for 
the corrections according to eq. ffTOj) with gi = 14.8(2) and guj = 2.9(2). The resulting data 
collapse is very good and only at very low temperatures corrections to scaling are stronger 



10 




FIG. 3: Scaling function ■d^ of the Casimir force in the 3D Ising model with (+— ) boundar 

conditions, compared with the mean-field prediction (solid line), the experimental data of ref. 
and the exact result for the two-dimensional Ising model (dashed line). 



try 

i 



and not fully accounted for by the ansatz in eq. ffTOj) . Such stronger corrections might be 
related to thejpresence of an interface in the system. "^-^ is compared with the experimental 



data in ref. 



9[, the prediction of mean-field theory lO] (solid line, normalized such that 



(MFT) , 



(0) = 'i9^^''(0)) and the corresponding result for the two-dimensional Ising model 



(dashed line). From the data set with L = 13 we estimate (0) = 5.97(2), in agreement 

with the experimental value i^^^^^iO) = 6(2) 9] but which is larger compared to the previous 
MC estimate ^9+_(0) = 4.900(64) (lOj and the analytical estimates ^^+P{0) = 3.16, 4.78. The 
latter depend on the approximant used to resum the field-theoretical e = 4 — rf-expansion up 
to 0(e) (see ref. 10|] for details). A good data collapse is obtained also for gi = 5.0(1) which, 

compared to fig. [3l yields an overall reduction of the amplitude of i)-^ by a factor R ~ 0.76. 

In addition, we expect the experimental data in ref. 9| to be affected by corrections to 
scaling already for x > 2, due to the relatively small corresponding value of < 30, where 
i ~ 3 A is the molecular scale in the experiment. In view of these difficulties, the agreement 
between the MC and the experimental data in fig. [3] is encouraging. In fig. H] the scaling 
function "(9++ for (++) BC has been obtained accounting for the corrections in eq. ( ITOl) with 
gi = 14.2(7) and g^ = 2.3(2). Using, instead, gi = 4.9(2) yields R ~ 0.77. Currently, for this 
BC and in film geometry no experimental data are available for comparison, but ^9++ can 



be compared with the prediction of mean-field theory 



101 ] (solid line, normalized as before) 
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FIG. 4: Scaling function "d+j^ of the Casimir force in the 3D Ising model with (++) boundary 
conditions, compared with the mean- field prediction (solid line) and the exact result for the two- 
dimensional Ising model (dashed line); • MC [lo| . 



and of the two-dimensional Ising model [12| (dashed line). From the data with L = 13 
we estimate 'i9++(0) = —0.884(16) which is slightly larger than the previous MC estimate 
'i9++(0) = —0.690(32) [10|] (indicated as a black dot in fig. HI still affected by finite-size 
corrections) and the corresponding field-theoretical predictions '*?+'+^(0) = -0.652,-0.346 
depending on the approximant used to resum the 0(e) series (see ref. 10|] for details). 
For PBC and L > 10 corrections to scaling turn out to be negligible and the result- 
ing function d-p in fig. [5] is in very good agreement with its previous determination 
based on the computation of the lattice stress tensor. The slight discrepancies 
might be due to the uncertainty in the normalization factor which had to be used in ref. . 
This agreement provides additional support concerning the reliability of our approach. Fig- 
ure [5] shows also the comparison with the available field-theoretical predictions [7| (solid and 
dashed lines) for x > up to 0{eY The discrepancies can be traced back to higher-order 



mg sea. 
in ref. 



nonanalytic contributions ~ e^^^ 22[ |. 



V. CONCLUSIONS 



We have presented a novel general approach to determine the universal scaling functions 
d of Casimir forces via Monte Carlo simulations. The corresponding results for the three- 
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FIG. 5: Scaling function "dp of the Casimir force in the 3D Isine; model with PBC. For comparison 
we report the normalized data set with L = 30 from ref. 8| which corresponds to the largest 
lattice size investigated therein. The solid (dashed) curve corresponds to the [1,0] ([0,1]) Fade 
approximant of the analytical prediction in refs. 3, 7]. 



dimensional Ising and XY UCs compare favourably with previous experimental, numerical 
and analytical results. Our predictions for and in a film geometry might also 
contribute to the understanding of Casimir forces acting on colloidal particles. 
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